rm(list = ls(all = TRUE))
setwd("C:\\Users\\lenovo\\Desktop\\R")
data(birthwt, package = "MASS")
str(birthwt)
library(dplyr)
birthwt <- birthwt %>%
  mutate(low = factor(low, labels = c(">=2.5kg","<2.5kg")),
         race = factor(race, labels = c("white", "black","other")),
         smoke = factor(smoke, labels = c("no smoking", "smoking")),
         ht = factor(ht, labels = c("no hypertension", "hypertension")),
         ui = factor(ui, labels = c("no uterine irritability", "yes")))
str(birthwt)
summary(birthwt)
library(epiDisplay)
summ(birthwt)
#5.1 数值型变量的描述性统计分析
cont.vars <- dplyr::select(birthwt, age, lwt, bwt)
length(cont.vars$age)
mean(cont.vars$age)
sd(cont.vars$age)
#同时计算多个变量的指定统计量
sapply(cont.vars,sd)
#计算偏度与峰度
#install.packages("psych")
library(psych)
describe(cont.vars)
#计算某个分类变量各个类别下的统计量
aggregate(cont.vars, by = list(smoke = birthwt$smoke), mean)
aggregate(cont.vars, by = list(smoke = birthwt$smoke), sd)
aggregate(cont.vars, 
          by = list(smoke = birthwt$smoke, race = birthwt$race),
          mean)
tapply(birthwt$bwt, INDEX = birthwt$smoke, mean)
#epiDisplay包的summ,按照分类变量绘制的有序点图
summ(birthwt$bwt, by = birthwt$smoke)
#psych包里的函数describrBy也可以分组计算与函数describe相同的统计量
describeBy(cont.vars, birthwt$smoke)
#利用dplyr包里的函数group_by和summarise来计算分组统计量
library(dplyr)
birthwt %>%
  group_by(smoke) %>%
  summarise(Mean.bwt = mean(bwt), Sd.bwt = sd(bwt))

#5.2 分类变量的列联表和独立检验
#5.2.1 生成频数表和列联表
#1 一维频数表
mytable <- table(birthwt$low)
mytable
#转换为百分比形式
prop.table(mytable)
round(prop.table(mytable)*100, 1)#保留一位小数
#epiDsiplay
library(epiDisplay)
tab1(birthwt$low)
tab1(birthwt$low, graph = "FALSE")#不输出条形图
#tab1函数可用于数值型变量来探索数据的缺失值与异常值
tab1(birthwt$age)
#2 二维列联表
mytable <- table(birthwt$smoke, birthwt$low)
mytable
#生成边际频数
addmargins(mytable)
#生成频率表
prop.table(mytable, margin = 1)#按行求比例
prop.table(mytable, margin = 2)#按列求比例
#用epiDisplay包里的函数tabpct
library(epiDisplay)
tabpct(birthwt$smoke,birthwt$low)
#3 多维列联表
mytable <- table(birthwt$smoke, birthwt$low, birthwt$race)
margin.table(mytable, 3)
margin.table(mytable, c(1,3))
addmargins(mytable)
#紧凑的格式输出
ftable(mytable)

#5.2.2 独立性检验
#卡方检验
mytable <- table(birthwt$smoke, birthwt$low)
mytable
chisq.test(mytable)
#Pearson's Chi-squared test with Yates' continuity correction
#对于每个格子的期望频数都大于5的样本，无需矫正
#期望频数表：
chisq.test(mytable)$expected
chisq.test(mytable, correct = FALSE)
#无论是否矫正，母亲吸烟情况与新生儿低体重都存在显著的关联

#Fisher精确概率检验
#一般用于总记录数n小于40、频数表中的某个期望频数很小
fisher.test(mytable)

#相对危险度与优势比
library(epiDisplay)
#第一个参数需要设置为结果变量
#cs求RR
#不知道为什么一直没有正常结果显示......
cs(birthwt$low, birthwt$smoke)
#cc求OR
cc(birthwt$low, birthwt$smoke)
mytable <- table(birthwt$low, birthwt$smoke)
cc(cctable = mytable)

#Cochran-Mantel-Haenszel卡方检验
#检验两个分类变量在调整（控制）第三个变量的情况下是否独立
#零假设：两个分类变量在第三个变量的每一层都是条件独立的
mytable <- table(birthwt$low, birthwt$smoke, birthwt$race)
mytable
mantelhaen.test(mytable)
library(epiDisplay)
mhor(mhtable = mytable)

#配对列联表的卡方检验
my.matrix <- matrix(c(11,2,12,3), nrow = 2)
mcnemar.test(my.matrix)

#5.3 连续型变量组间差异的比较 
#5.3.1 独立样本的t检验
#检查方差齐性
var.test(bwt~smoke, data = birthwt)
#t检验,默认方差不齐（var.equal = FALSE）
t.test(bwt~smoke, var.equal = TRUE, data = birthwt)
#另一种调用方法
group1 <- birthwt$bwt[birthwt$smoke == "no smoking"]
group2 <- birthwt$bwt[birthwt$smoke == "smoking"]
t.test(group1, group2, var.equal = TRUE)

#5.3.2 非独立样本的t检验
x <- c(0.84,0.59,0.67,0.63,0.69,0.98,0.75,0.73,1.20,0.87)
y <- c(0.58,0.51,0.50,0.32,0.34,0.52,0.45,0.51,1.00,0.51)
t.test(x,y, paired = TRUE)

#5.3.3 单因素方差分析
#Shapiro-Wilk正态性检验
tapply(birthwt$bwt, birthwt$race, shapiro.test)
#Bartlett方差齐性检验——对正态性较为敏感
bartlett.test(bwt~race, data = birthwt)
#Levene方差齐性检验——一种非参数方法
library(car)
leveneTest(bwt~race, data = birthwt)
#方差分析
race.aov = aov(bwt~race, data = birthwt)
summary(race.aov)
#Tukey法两两比较
TukeyHSD(race.aov)
#作图
plot(TukeyHSD(race.aov), las = 1)
#Bonferroni法、Holm法两两比较
#方法包括：
#“holm”, “hochberg”, “hommel”,
#“bonferroni”, “BH”, “BY”, “fdr”, “none”
pairwise.t.test(birthwt$bwt, birthwt$race, 
                p.adjust.method = "bonferroni")
pairwise.t.test(birthwt$bwt, birthwt$race, 
                p.adjust.method = "holm")


#5.3.4 组间差异的非参数检验
#明显偏态分布或者不具有方差齐性
wilcox.test(bwt~smoke, data = birthwt)

x <- c(0.84,0.59,0.67,0.63,0.69,0.98,0.75,0.73,1.20,0.87)
y <- c(0.58,0.51,0.50,0.32,0.34,0.52,0.45,0.51,1.00,0.51)
wilcox.test(x, y, paired = TRUE)

#各组之间相互独立，使用Kruskal-Waliis检验
kruskal.test(bwt~race, data = birthwt)
#各组之间不独立，如，重复测量设计，使用Friedman M检验
friedman.test()

#在控制第一类错误(拒真错误)的前提下
#使用Wilcoxon秩和检验进行多组间比较
#install.packages("PMCMRplus")
library(PMCMRplus)
comp <- bwsAllPairsTest(bwt~race, data = birthwt)
summary(comp)
plot(comp)

#5.4 用函数tablestack（）汇总双变量分析结果
str(birthwt)
attr(birthwt, "var.labels") <- c("Low birth weight", 
                                 "Mother's age(yr)", 
                                 "Mother's weight(lbs)",
                                 "Mother's race", 
                                 "Smoking status",
                                 "Number of premature births",
                                 "History of hypertension",
                                 "Uterine irritability",
                                 "Number of physician visits",
                                 "Birth weight(g)")
des(birthwt)
#low为结果变量；age、lwt、race、smoke为解释变量
tableStack(vars = age:smoke, by = low, dataFrame = birthwt)
tableStack(vars = 2:5, by = low, dataFrame = birthwt)
#用iqr指定检验方法（mother‘s wieght 改用 t test）
tableStack(vars = 2:5, by = low, iqr = lwt, dataFrame = birthwt)
#指定不同的输出形式，不显示检验的结果或检验的名称
tableStack(vars = 2:5, by = low, test = FALSE, dataFrame = birthwt)
tableStack(vars = 2:5, by = low, name.test = FALSE, dataFrame = birthwt)
#by可以是包含多水平的因子变量
tableStack(vars = c(low, age, lwt, smoke), by = race, dataFrame = birthwt)
#只显示total列
tableStack(vars = 2:5, by = "none", dataFrame = birthwt)
#用函数write.csv将函数tableStack的结果导出到.csv文件中
table1 <- tableStack(vars = 2:5, by = low, dataFrame = birthwt)
write.csv(table1, file = "table1.csv")
getwd()#显示当前的工作目录

#5.5 变量间的相关性
#5.5.1 连续型变量间的相关性
#1 Peason\Spearman\Kendall's Tau相关系数
library(dplyr)
cont.vars <- dplyr::select(birthwt, age, lwt, bwt)
cov(cont.vars)
cor(cont.vars, method = "pearson")
cor(cont.vars, method = "spearman")
cor(cont.vars, method = "kendall")

#2 相关系数的假设检验

#零假设：两变量之间不相关
#一个相关系数
cor.test(birthwt$lwt, birthwt$bwt)
#相关系数矩阵
library(psych)
#method默认为pearson，可用spearman、kendall
corr.test(cont.vars)
#各个相关系数的置信区间
print(corr.test(cont.vars), short = FALSE)

#3 偏相关
#在控制一个或多个变量时，另外两个变量之间的相关性
#install.packages("ggm")
library(ggm)
names(cont.vars)
#控制1等变量，看2、3
pcor(c(2, 3, 1), cov(cont.vars))
nrow(cont.vars)
r <- pcor(c(2, 3, 1), cov(cont.vars))
#偏相关系数的显著性检验
#q为条件变量的个数（要控制的变量的个数），n为样本量
pcor.test(r, q = 1, n = 189)

#5.5.2 分类变量间的相关性
library(vcd)
mytable <- table(Arthritis$Treatment, Arthritis$Improved)
#Phi系数只适合四格表
assocstats(mytable)

#配对列联表的Kappa统计量
my.matrix <- matrix(c(11, 2, 12, 33), nrow = 2)
library(epiDisplay)
kap(my.matrix)
#期望频数
chisq.test(my.matrix)$expected


#5.5.3 相关性的可视化
#数值变量
#散点图矩阵
pairs(cont.vars)
library(car)
scatterplotMatrix(cont.vars)

#相关系数矩阵
#install.packages("corrplot")
library(corrplot)
#圆面积的大小与颜色深浅代表相关性大小，蓝色正相关，红色负相关
corrplot(cor(cont.vars), tl.srt = 0)
#install.packages("corrgram")
library(corrgram)
corrgram(cont.vars, upper.panel = panel.pie)

#分类变量——关联图、马赛克图
library(vcd)
mytable <- table(Arthritis$Treatment, Arthritis$Improved)
mytable
#只能知道零假设是否被拒绝
chisq.test(mytable)
#关联图——展示二维列联表数据的一种工具
#矩形的高度与该单元格观测频数和期望频数之差成比例
#矩形的宽度与期望频数的成比例
assocplot(mytable)
#马赛克图——展示多维列联表
library(vcd)
mosaic(~ Sex + Treatment + Improved, data = Arthritis)


#同时展现数值型变量和分类变量的关联
#install.packages("GGally")
library(GGally)
dat <- dplyr::select(birthwt, age, lwt, bwt, race ,smoke)
ggpairs(dat)





